Decay of currents for strong interactions 
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The decay of current autocorrelation functions is investigated for quantum systems featuring 
strong 'interactions'. Here, the term interaction refers to that part of the Hamiltonian causing the 
(major) decay of the current. On the time scale before the (first) zero-crossing of the current, its 
relaxation is shown to be well described by a suitable perturbation theory in the lowest orders of the 
interaction strength, even and especially if interactions are strong. In this description the relaxation 
is found to be rather close to a Gaussian decay and the resulting diffusion coefficient approximately 
scales with the inverse interaction strength. These findings are also confirmed by numerical results 
from exact diagonalization for several one-dimensional transport models including spin transport in 
the Heisenberg chain w.r.t. different spin quantum numbers, anisotropy, next-to-nearest-neighbor 
interaction, or alternating magnetic field; energy transport in the Ising chain with tilted magnetic 
field; and transport of excitations in a randomly coupled modular quantum system. The impact of 
these results for weak interactions is finally discussed. 
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I. INTRODUCTION 

Perturbation theory is one of the main approaches to 
many-particle physics with a wide range of apphcations 
in the context of quantum transport, ranging from the 
investigation of Green's functions using Feynman graphs, 
the setup of a (quasi-)particle description by means of a 
Boltzmann equation, the derivation of the Green-Kubo 
formula within linear response theorji^, to many other 
applications in this context^. A particular application is 
the use of different projection operator techniquea^^^ for 
the realization of steady-state bath scenariosi2ri£ or the 
analysis of current autocorrelations^^-—^'. A perturbation 
theory is commonly employed due to the availability of a 
small parameter, e.g., weak particle-particle interactions, 
external scattering centers, or system-bath coupling. But 
additional assumptions are often required, such as the 
random phase and Markov approximatioiiSi'^. 
The concrete choice of a projection operator technique is 
a subtle task, whenever the addressed dynamics becomes 
non-Markovian and features memory effects^, typically 
occurring at short time scales. Because the relevant time 
scales are short in the case of strong perturbations, such 
non-Markovian effects appear in an already difficult case 
for any perturbation theory. However, the decay of the 
spin-current in the anisotropic Heisenberg chain^^'^^ at 
high temperatures has been well described for the case of 
large anisotropy parameters in Ref. [l^ by a lowest order 
prediction in the anisotropy, as obtained from a certain 
variant of projection operator techniques. Moreover, the 
resulting quantitative values for the diffusion coefficient 
have been brought into good agreement with numerical 
findings in the literatureiiii^ii^"— . But the perturbation 
theory in Ref. [l^ has focused on a single quantum model 
so far and has been carried out numerically on the basis 
of finite systems solely, leaving the origin of the observed 
agreements and disagreements as an open issue. Hence, 
one main intention of this paper is the extension of the 



perturbation theory to a wider class of quantum models 
and the analytical treatment of strong perturbations in 
the thermodynamic limit (and weak perturbations close 
to that limit). Furthermore, criteria for the validity of 
the lowest order prediction will be formulated and higher 
order corrections will be taken into account. By the use 
of these criteria and the comparison with numerically 
exact diagonalization (ED) the relaxation of the current 
is found to be well described in the lowest orders of the 
perturbation, even and especially if the perturbation is 
not weak. In particular the relaxation is rather close to 
a Gaussian decay and the resulting diffusion coefficient 
roughly scales with the inverse perturbation strength. 

This paper is structured as follows: In the next Sec. |TT] 
the general definition of the current and the connection 
between its autocorrelation and the diffusion coefficient 
is briefly reviewed at first. Then the perturbation theory 
for the decay of the current autocorrelation is introduced 
in Sec. Illll and the validity of the lowest order truncation 
for strong perturbations is discussed in detail here. In the 
following Secs. lIVIIVll the introduced perturbation theory 
is applied to several one-dimensional transport models 
in the limit of high temperatures, namely, the transport 
of excitations in a randomly coupled modular quantum 
system (Sec. IIVP : spin transport in the Heisenberg chain 
w.r.t. anisotropy, next-to-nearest neighbor interactions, 
different spin quantum numbers, or a staggered magnetic 
field (Sec. [V]); and energy transport in the Ising chain 
with a tiffed magnetic field (Sec. EH). The last Sec. IVHl 
closes with a summary and conclusion. 



II. CURRENT AND DIFFUSION COEFFICIENT 

In the present paper several (quasi-)one-dimensional 
and translationally invariant quantum systems will be 
studied, described by a respective Hamiltonian H . For 
such systems a globally conserved transport quantity X 
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will be considered, i.e., [H,X] = 0. The latter quantity, 
and the Hamiltonian as well, are both decomposable into 
N local portions Xr and hr-, corresponding to different 
spatial positions r: 
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Here, the Xr may be defined either exactly on the position 
of the hr , in between, or both. The above decomposition 
is further done in such a way that Heisenberg's equation 
of motion is of the form (ft = 1) 



dt 



— t[H, Xr] = l[hr- , Xr] + l[hr+ , Xr] 
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(2) 



where hr- and hr+ are located directly on the l.h.s. and 
r.h.s. of Xr, respectively. Because only the contributions 
from next neighbors r~ and r"*" are involved, this form 
may require the choice of a proper elementary cell. For 
instance, if additional contributions from next-to-nearest 
neighbors occur, a larger cell consisting of two or even 
more sites may be chosen. However, once a description 
in terms of Eqs. (dJ and ([2]) has been established, the 
local current is consistently defined by jr = i[xr,hr+^ 
and the total current reads 
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(3) 



r=l 



This paper will focus on the current autocorrelation 
function C{t) = (J(t)J(O)). Here, the time arguments 
of operators have to be understood w.r.t. the Heisenberg 
picture and the angles denote the equilibrium average at 
infinite temperature, i.e., essentially the trace operation: 
{...)— Tr{. . . l/dim'H. Particularly, the time-integral 
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(X2 



dt' C(i') 



(4) 



will be of interest. Apparently, for i — )• oo, the quantity 
X'(t) coincides with the diffusion constant according to 
linear response theory^'^. However, for any finite time, 
this quantity is also connected to the actual expectation 
value of local densities 



drit)=Tr{Xrit)piO)}-{Xr), 



(5) 



where p{0) represents an initial density matrix, featuring 
an inhomogeneous nonequilibrium density profile at the 
beginning. Concretely, the above quantity I?(t) and the 
spatial variance 



N 



Va.Tr{t) ^^dr{t)', 
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are connected via the relation^ 
d 



N 



dt 



Var,.(t) ^2V{t) 



(6) 
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Thus, whenever X'(t) is constant at a certain time scale, 
Varr(t) increases linearly at that scale, as expected for 
the case of diffusive dynamics. Contrary, — yields 
no increase (insulating behavior) and T^it) oc t leads to a 
quadratic increase (ballistic behavior). 
Strictly speaking, the relation in Eq. ([7]) is only fulfilled 
exactly for a class of initial states p(O)^'', representing 
however an ensemble average w.r.t. typicality'^®"'^*' or, 
more accurately, the dynamical typicality of quantum 
expectation values^. Hence, the overwhelming majority 
of all possible initial states p{0) is nevertheless expected 
to yield roughly a spatial variance corresponding to the 
ensemble average, if only the dimension of the relevant 
Hilbert space is sufficiently large. The latter largeness is 
certainly satisfied for most practical purposes. In other 
words, a concrete initial state p{0) is expected to fulfill 
approximately the relation in Eq. ([T]), at least if p(0) is 
not constructed explicitly to violate this relation. Note 
that the relation in Eq. (O remains still reasonable for 
lower temperatures, if all trace operations are performed 
including the statistical operator at a given temperature, 
see Ref. [23 for details. 

Generally, the dynamical behavior crucially depends 
on the considered time scale. At sufficiently short times 
transport generically is ballistic: The time-integral I?(i) 
firstly increases linearly, because the underlying current 
autocorrelation function C{t) has not decayed yet, at 
least not significantly, see Fig.[T] But, if the current is not 
strictly conserved, its autocorrelation function C{t) can 
further decay and may eventually reach a value close to 
zero at intermediate times. If C(t) additionally remains 
at a value in the vicinity of zero, its time-integral 'D{t) 
approximately becomes constant and finally develops a 
plateau at such time scales, see Fig. [TJ Then transport is 
diffusive and the respective diffusion coefficient is given 
by the height of the latter plateau. One main intention 
of this paper is to deduce the quantitative value of this 
diffusion coefficient perturbatively, especially in the case 
of strong perturbations, see the following Sec. IIIII 
However, even if I?(t) features a plateau at intermediate 
times, the dynamics does not necessarily stay diffusive 
for arbitrary long times. First, if the current is partially 
conserved, its autocorrelation function C(t) contains a 
non-decaying contribution (finite Drude weight). Such a 
contribution, independent of its weight, eventually leads 
to a renewed linear increase of I?(i) for sufficiently long 
times and consequently yields ballistic dynamics again, 
cf. Fig. [T] (red, dashed curve). Particularly, 'D{t) diverges 
for t ^ oo. The latter similarly happens, whenever the 
current autocorrelation function C{t) is (quasi-)periodic 
in t (finite recurrence time) . But the renewed increase of 
'D{t) in Fig. [T] (red, dashed curve) typically turns out to 
be a finite size effect, see Sec. IVland Refs. [l9l andls^ 
Second, apart from transitions towards ballistic behavior 
in the limit of long times, transitions towards insulating 
behavior can appear, of course. In that case the current 
autocorrelation function C{t) crosses zero and takes on 
a negative value, yielding a decrease of its time-integral 
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FIG. 1. (color online) Sketch for the possible time-dependence 
of (a) the current autocorrelation function C{t) and (b) its 
time- integral 'D{t). In (a) the current autocorrelation function 
C{t) may decay completely at short times and then remain 
zero at intermediate times (black, solid curve). At long times 
C{t) may decrease further to negative values (blue, dotted 
curve) or instead also increase to a positive value, due to a 
finite Drude weight (red, dashed curve). Such a Drude weight 
can already become dominant before the zero-crossing (green, 
dashed-dotted curve). In (b) the resulting time-integral T){t) 
may develop a plateau at intermediate times (black, solid 
curve), indicating diffusive dynamics. In that case the height 
of the latter plateau corresponds to the quantitative value 
of the diffusion constant V. The width of this plateau, and 
its existence as such, both depend on the concrete form of 
the underlying current autocorrelation function in (a) (other 
curves) . 



2?(i), see Fig. [T] (blue, dotted curve). However, all those 
transitions, occurring after a zero-crossing of the current, 
are not in the focus of this paper. In fact, the introduced 
perturbation theory in the following Section addresses 
the dynamics on the time scale before any occurrence of 
a zero-crossing. Strictly speaking, the theory itself does 
not even allow for a definite conclusion on the time scale 
after the first zero-crossing. But the theory yields at least 
an educated guess for I?(i), if C(i) is assumed to be more 
or less 'well-behaved', i.e., similar to Fig. [1] (black, solid 
curve). For instance, minor unsystematic fluctuations of 
the current autocorrelation function C (t) around zero do 
not contribute significantly to its time-integral I'(i) any 
further. Therefore, using this assumption, the diffusion 
constant is mainly set by the time scale before the first 
zero-crossing. For most of the concrete transport models 
in Sees. IIVIIVII this assumption is indeed supported by 
the numerical results from ED. 



III. LOWEST ORDERS TRUNCATIONS FOR 
STRONG PERTURBATIONS 

In this Section the perturbation theory for the analysis 
of the dynamical behavior of the current autocorrelation 
function C{t) is introduced. However, this theory is not 
restricted exclusively to the autocorrelation function of 
the current and may be applied analogously to any other 
observable of interest. Moreover, even the decomposition 
of the Hamiltonian according to Eq. ([IJ, considered for 
the definition of the current only, is not necessary in the 
following. 

Generally, a strategy for the perturbative description of 
the dynamical behavior of an autocorrelation function, 
such as C{t), is the application of projection operator 
techniques, particularly, the time-convolutionless (TCL) 
method^iiSiii^, as used here. This particular method, and 
the well-known Nakazima-Zwanzig (NZ) method^"^ as 
well, are both applied commonly in the context of open 
quantum systems. In the context of closed quantum sys- 
tems the Mori-Zwanzig memory matrix formalis m^i^i^^i^^ 
is applied more commonly. The latter formalism is, how- 
ever, similar to the NZ method. 

For the application of projection operator techniques a 
suitable projection (super-)operator V has to be defined 
firstly, projecting a density matrix p(t) onto a relevant 
subspace. This relevant subspace has to include at least 
the identity 1 and the observable of interest J. Thus, a 
natural choice of V is given by 



Vpit) 



1 



dimH 



{J Pit)) 
(J2) 



J 



(8) 



Here, the Ai denote additional observables which are of 
interest by themselves or crucially affect the dynamical 
behavior of the single observable of interest J. Without 
loss of generality, the set of all operators {l,J,Ai) may be 
assumed to be orthogonal w.r.t. the trace operation, i.e., 
(J) = {A.,) = {JAi) = (AiAj) =0. If J represents the 
current, 1 and J are always orthogonal, since no current 
flows in equilibrium, i.e., (J) = 0. This orthogonality and 
the normalization factors in Eq. ([5]) guarantee — V, 
i.e., V satisfies the property of a projection. Moreover, 
for initial conditions p(0) in the span of I and J, the 
two properties (I — V)p{0) — and {JV p{t)) oc C{t) are 
fulfilled. Especially the latter property demonstrates the 
connection to the autocorrelation function C (t) . 

Apart from the definition of an appropriate projection 
(super-)operator, the application of projection operator 
techniques requires the decomposition of the Hamiltonian 
into the form H = Hq + AV with Hq as the uncoupled 
system and with V as the 'interaction'. The parameter 
A is the strength of the interaction. This decomposition 
is independent from the one in Eq. ([1]), considered for 
the definition of the current solely. The decomposition 
of H is usually done in such a way that the interaction 
strength A is a small parameter and the observables of 
interest are strictly conserved in the uncoupled system. 
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e.g., [Hq, J] = Oiiii^. However, if the interaction strength 
A is a large parameter, it will be sufficient to assume that 
the dynamics in Hq is much slower than in A T^. For very 
strong interactions also the eigensystem of Hq will not be 
necessary. But for weak interactions the eigensystem is 
indispensable and has to be found analytically or at least 
numerically. 

Once the projection has been chosen and the uncoupled 
system has been identified, the TCL formalism routinely 
yields a closed and time-local differential equation for the 
evolution of V pi (t) in the interaction picture^ 

^ r pi{t) - g{t) V pi{t) + i{t) (1 - V) p{o) (9) 

and avoids the often troublesome time-convolution which 
appears, e.g., in the context of the NZ method. For initial 
conditions p{0) with (1 — V) p{Q) = the inhomogeneity 
on the r.h.s. of Eq. ^ vanishes, e.g., for the considered 
p(0) in the span of 1 and J. The generator Q{t) is given 
as a systematic perturbation expansion in powers of the 
interaction strength A, namely, 

oo 

g(0 =^A'g,(t), fe-iW = 0. (10) 

i=l 

The odd contributions of this expansion vanish in many 
situations and also for all concrete transport models in 
the subsequent Sections. Hence, the lowest nonvanishing 
contribution is the second order G2it). It reads, using the 
Liouville operator C{t) = —i[Vi{t),pi{t)], 

Q2{t) = / AhVC{t)C{t^)V . (11) 
Jo 

The next lowest nonvanishing contribution is the fourth 
order Qiit). It is given by 

nt I'll nt2 

Gi{t)^ / dti / dt2 / dt3 
Jo Jo Jo 

VC{t)C{h){l-V)C{h)C{h)V 
-V L{t) C{t2)V L{h) C{h)V 
^VC{t)C{h)VC{h)C{t2)V. (12) 

Generally, Eq. ^ leads to a rate matrix equation for 
the evolution of the expectation values {JV pi{t)) and 
{Ai V pi{t)) of all operators in the chosen projection. For 
the remainder, however, the discussion will focus on the 
projection onto a single operator J. In that case Eq. © 
yields a rate equation for the dynamical behavior of the 
expectation value {JV pi{t)) cx {J{t) Ji{t)) only. As long 
as the time evolution of J w.r.t. Hq is negligibly slow, 
say, Ji{t) ~ J, the latter expectation value is identical to 
the autocorrelation function C{t). Hence, Eq. (|9]) can be 
rewritten as 

Cit) = -[A2 R,{t) + A^ R,{t) + ...] C{t) (13) 



with scalar rates A^i?2(0 and A'*i?4(t), resulting from 
Eqs. (|10|) -(fT2 |) . Concretely, the second order rate R2{t) 
is given by 

^{d,. /(,,), /(M^™™ ,14) 

and the fourth order rate i?4(t) reads 

Riit) = dh dt2 / dts 
Jo Jo Jo 

f{t-ti) f{t2~h) 
+ f{t-t2) I{tl~h) 

+f{t~h) f{h-t2) 
{Uv,{t,)],Vi{t)]Uv,{t3)],Vi{t2)]} 

m ■ ^ ^ 

The above Eqs. (|13|) - (fT5|) build the framework for the 
analysis of the decay of the current autocorrelation C{t) 
in the next Sections. 

So far, the rate equation is formally exact, at least if 
the rates in all orders of the perturbation expansion are 
taken into account. But already the concrete evaluation 
of the fourth order rate is typically a highly nontrivial 
task, both analytically and numerically. The rest of this 
Section will therefore discuss the possibility of a second 
order truncation in detail. The quality of this truncation 
may crucially depend on the choice of the projection and 
is commonly expected to be justified only in the limit of 
weak interactions^. 

In the limit of weak interactions the truncation to lowest 
order usually relies on the fact that decay times become 
arbitrary long, if only the strength of the interaction is 
sufficiently small. In the case of long time scales, e.g., in 
the Markovian limit the rates in all orders are assumed to 
take on constant values A* Ri. Due to the independence 
of the rates from time, the contribution of the ith order 
rate is essentially determined by the overall scaling factor 
A* solely. Therefore, in the limit A — > 0, the dominant 
contribution is given by the lowest order rate A^ R2. By 
the use of this assumption, Eq. ([T5)) directly predicts the 
purely exponential decay C{t) = exp(— A^i?2 1) ( J^). For 
times above the decay time r — l/(A^i?2) the diffusion 
coefficient according to Eq. ^ consequently becomes 

1 1 (J2) 

25wcak(i) = -^^ = const. , (16) 

i.e., the diffusion coefficient scales as 1/A^, if the other 
quantities (J^) and {X"^) do not depend on A. However, 
the underlying assumption of constant rates may not be 
fulfilled, or at least not for the projection onto a single 
observable. A possible criterion for the validity of such a 
lowest order prediction might be the independence from 
the number of observables, see Sec. |Vl 
In the limit of strong interactions a truncation to lowest 
order appears to be not reasonable. On the one hand the 
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interaction strength is large as such, and on the other 
hand the relevant time scales are short and the rates in 
all orders are explicitly time-dependent here, e.g., in the 
non-Markovian limit. Moreover, in that limit significant 
differences between the lowest order predictions of the 
TCL and NZ approach are commonly expected due to 
memory effects^. But, as a consequence of the explicit 
time-dependence of rates, a truncation to lowest order 
will turn out to be justified anyhow. 
To this end consider the rates in Eq. (ITil) and (|15l) . For 
sufficiently short times the integrands can be treated as 
approximately time-independent and the time-integrals 
can be easily performed. Thus, in the case of short time 
scales, Eq. ((Tl)) and ((T5|) are approximated by 

R2{t)^r,t,r,^^-^1:^, (17) 

i?4(i) ^rj^,u = 3rl- ^^^'^'^^p'^ , (18) 

where the time-dependence of the ith order rate appears 
only as the power Apparently, the approximations 

do not contain any feature of the uncoupled system Hq, 
neither its eigenvectors nor its eigenvalues. But Hq still 
determines the maximum time for the validity of these 
approximations. Or, in other words, the approximations 
are valid, if the resulting decay time r becomes shorter 
than typical decay times in Hq, see below. 
Now assume for the moment that a truncation to lowest 
order was indeed correct. Then Eq. (13) directly predicts 
the strictly Gaussian decay C(t) = exp(— r2 t^/2) ( J^) 
and for times above the decay time r = -^2/ (A y^) the 
diffusion coefficient in Eq. thus becomes 



2?strong(i) = ^ Y ^ J]^ = const. , (19) 

i.e., the diffusion constant scales as 1/A in that case, in 
contrast to the former case of weak interactions. Such a 
prediction is only correct, if higher order rates represent 
minor corrections at the relevant time scale, i.e., up to 
the decay time r. Hence, the ratio of the fourth to the 
second order rate may be compared at this particular 
point in time. According to Eqs. (IT7)) and ([T5)) . the ratio 
reads 

A2 i?2(T) 6r2 3(i[J,F]2)2 ^ 

and does not depend on the interaction strength A. The 
latter independence from A also holds true for ratios with 
higher order rates. Thus, the validity of the second order 
truncation for strong interactions is connected to static 
expectation values, involving exclusively the form of the 
interaction and not its strength as such. In a sense the 
situation is rather similar to weak interactions. 
In principle the ratio in Eq. (j20p can certainly be a large 
negative number. But in many situations and also for all 



concrete transport models in the following Sections the 
ratio turns out to be rather close to zero. Due to such a 
small ratio, the relaxation of the current is found to be 
well described in terms of the second order prediction, at 
least up to the decay time. This restriction is necessary, 
simply since the introduced approach does not allow to 
describe zero-crossings of the current, see Eq. (jl3p . Such 
a zero-crossing essentially requires the divergence of the 
whole series expansion at this point in time. However, it 
will be demonstrated that the relaxation of the current is 
described almost exactly up to the (first) zero-crossing for 
strong interactions, if fourth order corrections are taken 
into account. Particularly, the case of strong interactions 
will be shown to include physically relevant and not only 
pathologic situations. 



IV. MODULAR QUANTUM SYSTEM 

This Section will investigate the transport of a single 
excitation amongst local modules in an one-dimensional 
configuration, as a first example for a concrete quantum 
systemi^i2£ii2^. This modular quantum system is one of 
the few models which have been reliably shown to ex- 
hibit purely diffusive transport. In fact, pure diffusion 
has been found from different approaches and also the 
same quantitative value of the diffusion coefficient has 
been deduced. But all these approaches have focused on 
the case of weakly coupled modules. The present inves- 
tigation will address the so far unexplored case of strong 
couplings. 

To start with, each of the N local modules features an 
identical spectrum, consisting of n equidistant levels in a 
band with the width 6e. Therefore the local Hamiltonian 
at the position r is given by hr — J2f_i l^^^/n ^A 
the one-particle basis \r,^). The next- neighbor coupling 
between two local modules at the positions r and r + 1 
is supposed to be 

Av^ = A^c^,^|r,/i)(r + l,i^|+H.c. (21) 

with the overall coupling strength A. The r-independent 
coefficients c^_i/ are complex, random, and uncorrelated 
numbers: their real and imaginary part are both chosen 
corresponding to a Gaussian distribution with the mean 
and the variance 1/2. Here, a single realization of these 
coefficients is considered (and not an ensemble average 
over different realizations). The total Hamiltonian reads 
H = Hq -t- A 1/, where Hq denotes the sum of all local 
Hamiltonians hj. and V represents respectively the sum 
of all next-neighbor couplings Vr- 

Of particular interest is the probability for finding the 
excitation somewhere within the rth local module. Such 
probabilities correspond to local density operators of the 
form Xr — lA- Their total sum is the identity 

1 and a globally conserved quantity. Therefore, according 
to the scheme in Sec. [TTl the definition of the associated 
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local current is given by jV = ^[Pri ^Vr], yielding 

jr = Aj2i c^,^ \r,^l){r + 1,1^1+ H.c. (22) 

with a form similar to Vr in Eq. ([21]) . The total current 
J does not commutate with H, Hq, or V. 

Although the modular quantum system is not meant to 
represent a concrete physical situation, it may be viewed 
as a very simplified model for a chain of coupled atoms 
or molecules. In that case the hopping of the excitation 
from one module to another corresponds to transport of 
energy. Alternatively, the modular quantum system may 
be illustrated as an idealized model for non-interacting 
particles on a more-dimensional lattice. In that case the 
hopping of the excitation between modules corresponds 
to transport of particles between chains or layers. There 
also is a relationship to the three-dimensional Anderson 
model^i^^, even though the modular quantum system 
does not contain any disorder, despite the random choice 
of the coefficients c^^i, in Eq. ([2T|) . 

However, since the current J does not commutate with 
the uncoupled system Hq, the autocorrelation function 
C{t) can not be analyzed by the introduced approach in 
Sec. lIIIi if the coupling strength A is small. But for such 
small A the dynamical behavior of the autocorrelation 
function C (t) can be analyzed in a different way. To this 
end consider the total Hamiltonian H = Hq + A for the 
case of sufficiently small A. In that case the eigensystem 
of H is essentially given by Hq, i.e., the eigenvectors are 
|r, /x) and the eigenvalues are /iSe/n. Obviously, C{t) can 
be determined more or less exactly: Since the spectrum 
features the width Se, the autocorrelation function C{t) 
fully decays at a first time scale tq ~ n/Se and eventually 
recurs completely at a second time scale Tq = 2'Kn/de, 
simply due to the equidistant levels of the spectrum. But, 
within the possibly wide time window between these time 
scales, C{t) remains zero and thus the diffusion constant 
I?(i) according to Eq. (jH) becomes constant. Concretely, 
the diffusion constant reads-^ 

_ , , TT (J2) 27r A^n 

2?woak i = Y- = — I = const. 23 

for sufficiently many levels in the spectrum. The scaling 
factor A^ results from ( J^) = 2A^n. This result coincides 
with Fermi's Golden Rule and has been obtained already 
from different approache a^^'^°i'^^ . 

In the so far unexplored case of strong couplings the 
approach in Sec. Illll becomes applicable. The application 
only requires that the decay due to A proceeds much 
faster than the decay due to Hq, i.e., C{t) has to decay 
at a time scale t < tq. If A is large, this requirement is 
naturally fulfilled. In fact, at such a time scale below tq 
no correlation function w.r.t. Hq has decayed yet. Hence, 
the second order rate in Eq. ([14]) can be well described 
by the linear approximation in Eq. (|17p , see Fig. [5| The 
use of {i[J,Vf) SA^n^ and (J^) = 2A2n concretely 
leads to r2 = 4 n. The resulting second order prediction 
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FIG. 2. (color online) The second order decay rate R2{t) for 
the excitation current in the modular quantum system (black, 
solid curve). Parameters: n = 100 and Se — 0.5. This decay 
rate is well described by the linear approximation Ant at a 
time scale below n/St (red, dashed line). 

for the diffusion constant in Eq. (|19l) consequently is 

^ , , 1 / TT LP) Inn ^ 

I^W0 = ^^3;^^ = ^-A = const. (24) 

above the relaxation time r = l/(\/2n A). In Fig. [3]the 
analytical second order prediction is compared with the 
direct numerical result for the diffusion coefficient by the 
use of ED. Since for the modular quantum system the 
linear growth of the Hilbert space can be compensated 
by the translation symmetry, rather many modules are 
treatable. Apparently, the agreement between the second 
order prediction and numerics is very good, despite the 
limit of strong interactions. The minor deviations can be 
further reduced by fourth order corrections in terms of 
the cubic approximation in Eq. ([TSj) . Concretely, the use 
of ([[J, V], VY) = 80 A2 leads to r4 = 8 n^. With such 
fourth order corrections the agreement in Fig. [Hj becomes 
excellent. Remarkably, the corrections are small, because 
the ratio in Eq. (|20p takes on the value 1/6. 



V. HEISENBERG CHAIN 

In this Section spin transport in the Heisenberg spin 
chain will be investigated, as another concrete example 
for an interacting man?/-particle quantum system, going 
beyond the simplified single-particle model in the last 
Section. This investigation will firstly focus on a certain 
generalization of the standard Heisenberg spin- 1/2 chain, 
taking into account the effect of anisotropic nearest and 
next-to-nearest neighbor interactions. The Hamiltonian 
concretely reads H — Hq + AV^, where the operators Hq 
and V are given by 

i^o = Jec^s?:<+i + .s^s:{+,, (25) 

r 

V^J,,Y.'>r+l + ^slsl^2- (26) 
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FIG. 3. (color online) The time-dependent diffusion constant 
T>{t) for the modular quantum system with strong interactions 
(a) A = 0.1, (b) A = 1.0. Parameters: N = 1000, n = 100, 
and 5e — 0.5. The theoretical predictions according to TCL2 
(dashed, red curves) are already in good agreement with the 
numerical results from ED (solid, black curves), while minor 
corrections are given by the approximation of TCL4 (dotted, 
blue curves). As predicted by theory, 2?(t) takes on a constant 
value 2? oc A at intermediate times above the relaxation time 

T CX 1/A. 



Here, Jcc is the exchange coupling constant, A refers to 
the anisotropy parameter, the matrix s* represents the 
ith component of the spin-1/2 operator at site r, and 
the parameter A2 specifies the strength of an additional 
next-to-nearest neighbor zz-interaction. For A2 = the 
Hamiltonian obviously reduces to the usual anisotropic 
Heisenberg spin-1/2 chain (XXZ model). 
The transport of spin or magnetization corresponds to 
local density operators Xr = s^. Their sum is and a 
globally conserved quantity. Therefore, according to the 
scheme in Sec. the associated current can be written 
in the well-knovifn form (see, e.g., the reviews [2TI and [23) 



J=Je 



(27) 



and commutates with Hq. In particular, (J^) = 
and {X'^) — N /A, where the trace operation is performed 
over the full 2^-dimensional Hilbert space. In fact, the 
following investigation will not be restricted to a specific 
M-subspace of . However, the dominant contribution 
to transport stems from the largest subspaces around 
M — {) ('half filling'). The dynamics in these subspaces 
can be diffusive, while the dynamics in the subspaces 
with \M\ ^ ('dilute filling') is expected to be ballistic, 
see below. 



Since the eigensystem of Hq is indispensable for the 
application of the introduced approach in Sec. Illli it is 
convenient to perform the Jordan- Wigner transformation 
onto spin-less fermions--- firstly and an additional Fourier 
transformation afterwards. The operators iJo, V ^ and J 
in Eqs. (I^5|) - (P7)) can then be rewritten as^ 

Hq = Jcc ^£knk , ek = cos{k) , 

k 



k,l,q 



dek 



J ^ Jcc^^ jk rik , jk = ^ = - sin(A:) 

k 



(28) 



Here, = a[,afc denotes the particle number operator 
for a spin-less fermion with the momentum k — 2ni/N 
and is written as the product of respective creation and 
annihilation operators. In this picture, Hq describes the 
dispersion of non-interacting particles, while V is the 
interaction between two particles, located at nearest or 
next-to-nearest sites, and J plays the role of a particle 
current. Since Hq and J are both diagonal, J is strictly 
preserved in the absence of V and also in the one-particle 
subspace. As a consequence a single particle propagates 
ballistically. 

If additional next-to- nearest neighbor xx- / yy-tevvas are 
added to Eq. (|25l) [and hence to Eq. ([ST)) ], such a picture 
can not be established: Hq does not become diagonal by 
the use of the Jordan- Wigner transformation and J does 
not commute with Hq . But these facts do not imply that 
the approach as such becomes not applicable. Hq can be 
diagonalized at least numerically and, for large A, the 
commutation of J with Hq is not required in the strict 
sense. Obviously, the situation is similar to the modular 
quantum system in Sec. |lVl However, it turns out that 
the decay of J w.r.t. Hq is comparatively fast. This fast 
decay restricts the applicability of the approach to very 
large A, i.e., close to the less interesting Ising limit. Thus, 
a situation with next-to-nearest neighbor xx- / yy-tevvas 
will not be discussed further. 

For operators of the form in Eq. ((28|) an exact analytical 
formula for the second order decay rate R2 (t) in Eq. 
can be derived. In fact, the derivation of such a formula 
only requires the concrete evaluation of the expectation 
value {[nk,V]i{t)[ni,V]) in the interaction picture, i.e., 
w.r.t. Hq. Even though the concrete evaluation appears 
to be a straightforward task at the first view, it turns out 
to be a both subtle and lengthy calculation. Nevertheless, 
after such a calculation the latter expectation value can 
be finally given as 

{i[nk,V]i{t)i[nuV]) ^ 

(J2) N'^ ^ 

Sk,l ^{Vk-q ~ Vm^qf C0s[(efc + e„i - f-k+rn-q - f-q) A 

■m 

+ {vk-q - Vl-qf COs[(efc + £; - €k+l-q - ^q) t] 



-2 (vk-q - vi-kf cos[{eq + ei- eq+i-k - efc) t] 



(29) 
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FIG. 4. (color online) The second order decay rate R2{t) for 
the spin current in the anisotropic Heisenberg chain without 
additional terms (black, solid curve) and with an additional 
next-to-nearest neighbor z 2- interaction of the same strength 
A2 = A (red, dashed curve). Parameters: N = 2000 [possible 
due to Eq. (|30p ]. These decay rates are well described by the 
linear approximation t/2 at a time scale below 1/Jcc (blue, 
dotted line). Finite size effects do not occur at time scales up 
to 500/ Jcc. 
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with {ife = Re Vk ■ Consequently, for a linear combination 
of the form Ai = ^Ylik "^fe "•fc '^^'^ finds 




including the second order decay rate i?2 (i) for i = j = \ 
and a\ = jk- The more general notation in Eq. pOI) will 
become useful later. Because this equation only involves 
a sum over three momenta k^l, and q, it can be evaluated 
numerically for several thousands of spins and, say, in 
the thermodynamic limit. Concretely, N = 2000 will be 
chosen in the following. For that choice finite size effects 
do not appear at time scales up to 500/ Joe- For instance, 
such finite size effects occur at time scales on the order of 
10/ Joe, if TV = 20 is chosei>i2, e.g., the maximum number 
of spins for ED. 



A. The case A2 = A 

First, the case of additional next-to- nearest neighbor 
zz-interactions of the same strength may be discussed 
in detail, i.e., A2 — A. This particular case appears to 
be less controversial, since Drude weights are commonly 
expected to vanish in the thermodynamic limit due to 
non-integrability^i"— , at least if A2 (= A) does not be- 
come too small. For the case A2 = A the second order 
decay rate R2 (t) in Fig. U indeed takes on a form, as 
already considered in Sec. IIIII i?2(i) firstly increases lin- 
early at short time scales below 1/Jcc and then becomes 



FIG. 5. (color online) The time-dependent diffusion constant 
T>{t) for spin transport in the anisotropic Heisenberg chain 
with an additional next-to-nearest neighbor zz-interaction of 
the same strength A2 — A for large anisotropy parameter 
(a) A — 1.0 and (b) A = 2.0. The theoretical predictions of 
TCL2 for A'^ = 2000 (red, dashed curves) are already in good 
agreement with the numerical results from ED for TV — 10, 
18, and 20 (black, solid curves). Additional corrections on the 
order of 20% are given by the approximation of TCL4 (blue, 
dotted curves). In (a) the approximation of TCL4 is a slight 
overestimation. 

constant at longer time scales. Thus, the second order 
predictions for the diffusion coefficient V can directly be 
formulated according to Eqs. ((TB| and (IT^ . By the use 
of R2 ~ 0.71 Joe from Fig. 3] and r2 — 1/2 J2^ these pre- 
dictions read 

I'woak ~ 0.70 , I?stiong = • (31) 

For anisotropy parameters above A ^ 1 the decay of the 
current autocorrelation function takes place at a short 
time scale, i.e., where i?2(i) scales more or less linearly 
with time. Since Drude weights for such A are already 
sufficiently small for TV ~ 20, a direct comparison with 
the numerical results from ED becomes possible here, as 
shown in Fig. [S] By the use of the exact R2{t) in Fig. 0] 
the second order predictions for Vit) are already in good 
agreement with ED. Additional corrections on the order 
of 20% are given by the fourth order approximation in 
Eq. p^ . For A = 1 the latter approximation for short 
times seems to be a slight overestimation and may not 
be used further for smaller A, i.e., when the relevant 
time scales for the decay of the current autocorrelation 
function become longer. 

Since Drude weights always become dominant for small 
A and finite systems, a direct comparison between the 
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FIG. 6. (color online) The time-dependent diffusion constant 
©(t) for spin transport in the anisotropic Heisenberg chain 
with an additional next-to-nearest neighbor zz-interaction of 
the same strength A2 = A = 0.5. The theoretical prediction 
of TCL2 for A'' = 2000 (red, dashed curve) is consistent with 
the numerical results from ED for — 10, 18, and 20 (black, 
solid curves). Another projection with more observables leads 
to a modified prediction TCL2+ for A' = 200 close to TCL2 
(blue, dotted curve). 
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t [units of 1/J 



20 



second order prediction and the numerical results from 
ED is difficult in that case. But for A = 0.5 theory and 
numerics are at least consistent, see Fig. |6l Because the 
fourth order approximation is not available, the validity 
of the second order prediction may be confirmed by its 
independence from the chosen projection. To this end the 
projection may be extended to the full diagonal space, 
consisting of linear combinations Ai of particle number 
operators Uk- As outlined in Sec. IIIIl such a extension of 
the projection yields a rate matrix equation, namely, 

^^{.mA,) =-A^J2^2'it){Jit)Aj) (32) 
j 

with decay rates i?2 "' (t) according to Eq. ([501 . This rate 
matrix equation can be solved numerically by standard 
algorithms for, e.g., N = 200. This modified prediction 
of, say, TCL2-I- for A = 0.5 turns out to be rather close 
to the original prediction of TCL2, see Fig. [S] Because 
both predictions become identical for larger A, TCL2-I- 
is not indicated explicitly in Fig. [5l However, TCL2 and 
TCL2-I- begin to differ significantly for smaller A and the 
validity of Eq. pip in the limit of very weak interactions 
is questionable without the consideration of higher order 
decay rates. 

B. The case A2 = 

The case A2 = is rather controversial due to the 
intcgrability of the Hamiltonian in terms of the Bethe 
Ansata^. In particular, for the isotropic point A = 1, 
there still is an unsettled debate about the finiteness of 
Drude weights in the thermodynamic limit: While Drude 
weights are widely expected to vanish for A > 1— i^, they 
may already become zero for A = 1—"— but definitely 



FIG. 7. (color online) The time-dependent diffusion constant 
©(t) for spin transport in the anisotropic Heisenberg chain 
for the anisotropy A = 1.5 and the spin quantum numbers 

(a) s = 1/2 and (b) s = 3/2. The numerical results from 
ED (black, solid curves) are shown for (a) A' — 10, 18, 20 and 

(b) A' = 7, 8, 9. The TCL2 predictions (red, dashed curves) 
are based on the data in Fig. U and Ref. [l9l . respectively. No 
TCL4 corrections (blue, dotted curve) are required in (b). 



non-zero for all < A < 1—, see also Ref. [3. How- 
ever, in the present approach the situation is found to 
be similar to the previous case A2 = A, see Fig. S) The 
second order decay rate i?2(i) is almost the same, i.e., 
with a slightly reduced value R2 ~ 0.58 Joe- Remarkably, 
the latter value can already be supposed on the basis of 
N ^ 20, i.e., Hq in Eq. ^ contains only 10 different 
energiesli^.. Due to Fig. IH the second order prediction 
for I?strong in Eq. pi]) remains unchanged and I?woak be- 
comes 



2?wcak ~ 0.86 . (33) 

But this prediction has to be considered carefully, since 
it depends on the chosen projection in the limit of very 
weak interactions, analogously to the case A2 = A. In 
fact, significant differences between TCL2 and TCL2-f' 
occur already for A ~ 1. Nevertheless, for larger A the 
TCL2 prediction is again found to be in good agreement 
with the numerical results from ED, see Fig. [7] (a). The 
additional incorporation of the TCL4 approximation in 
Eq. goes beyond Ref. [l^ and explains the reported 
deviations from ED and other approache s 1 1 ^ 1 ^3 ^ g^jg^ 
the perturbative approach in Ref. [l6l 
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(b) A = 2.0 (5B/J^^ = 0.25) 
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FIG. 8. (color online) The time-dependent diffusion constant 
T>{t) for spin transport in the anisotropic Heisenberg chain in 
the presence of an alternating magnetic field SB/Jcc = 0.25 
for (a) A — 1.0 and (b) A = 2.0. The numerical results from 
ED (black, solid curves) are shown for A'' — 10, 16, and 18. 
The theoretical predictions of TCL2 (red, dashed curves) are 
also based on numerical data for N = 18, even though TCL2 
does not depend on A'^ any further. 



C. The case s > 1/2 (A2 = 0) 

For spin quantum numbers s > 1/2 the operators Hq, 
V, and J in Eqs. p5|) - (l27l) are formally identical. But 
in that case Hq can not be brought into diagonal form 
by the use of the Jordan- Wigner transformation. The 
latter lack of a diagonal form is not a substantial draw- 
back, since J does not commutate with Hq for s > 1/2. 
Thus, the investigation is anyway restricted to the limit 
of strong interactions, similarly to the modular quantum 
system in Sec.lIVI Such an investigation has already been 
done in detail in Ref. [l^. For illustration, however, an 
example for s = 3/2 and A = 1.5 is shown in Fig. [7] 
(b). In a sense it is intriguing to see that the agreement 
between the theoretical prediction of TCL2 and the nu- 
merical result from ED is best, if [J, Hq] ^ 0, see also 
Fig. [31 It is worth to mention that the investigation in 
Ref. lis suggests that, at high temperatures, the diffu- 
sion constant scales with the spin quantum number as 
I'strong (X \/ s{s -|- 1), scc also Rcfs. [2J-[2^ This scaling 
supports classical simulations at high temperatures^i^. 

D. Alternating magnetic field (A2 = 0) 

So far, the investigation at high temperatures does not 
depend on the presence of a homogenous magnetic field 



B in z-direction, i.e., it does not change for an addi- 
tional Zeeman term Hb = BS'' in Eqs. (I25|) and (|26l) . 
But the situation changes in the presence of an inho- 
mogeneous magnetic field^^— , e.g., for an alternating 
sequenceiiii^ii^'^ 

H,.^^.±S^,i, ,34, 

r 

where SB denotes the strength of the alternation. Since 
[J, Hb] 7^ 0, a magnetic field of this form represents an 
additional scattering mechanism and has to be added to 
V in Eq. (I26p . In principle the alternating magnetic field 
in Eq. p4|) may also be written in the representation of 
spin-less fermions and an exact analytical formula for the 
second order decay rate R2{t) may again be derived. But, 
because Eq. p4|) obviously is no two-particle interaction, 
i?2(i) is not of the form in Eq. (PH)) . However, the case 
of an alternating magnetic field is primarily considered 
here in order to demonstrate potential difficulties of the 
approach at hand. For that reason the decay rate i?2(i) 
is directly evaluated numerically for iV ~ 18 by the use of 
EEiiiii^, e.g., finite systems of this size are sufficient for 
the limit of strong interactions^^. In Fig. |8]the resulting 
TCL2 prediction for 2?(t) is shown for SB/J^c = 0.25 and 
different A. Apparently, there still is a good agreement 
with the numerical results from ED apart from the, say, 
oscillation in Fig.[8](b). This oscillation takes place after 
a zero-crossing of the underlying current autocorrelation 
function^^ and is therefore not captured by the present 
approach. Because the amplitude of such oscillations is 
known to increase with the strength of the alternation^, 
the TCL2 approach yields only meaningful predictions 
for not too large 6B. 

VI. ISING CHAIN 

This Section will deal with another concrete quantum 
system which also allows to clarify potential difficulties of 
the approach at hand. This quantum system is an Ising 
spin-1/2 chain in the presence of a, say, tilted magnetic 
field. The Hamiltonian concretely reads H — Hq + B^V, 
where the operators Hq and V are given by^^ii^i^i^ 

Ho=B,S^ + J,,J2s',s',+,,V = S\ (35) 

r 

where Bi and are the ith component of the magnetic 
field B and total spin S, respectively. For instance, one 
might think of a magnetic field which was originally in 
line with the z-direction and has been rotated about the 
y-axis with the angle a = SLictaii^Bx/ Bz). 
Since [S^, H] for this model, spin or magnetization is 
not a suitable transport quantity here. However, energy 
is always an appropriate transport quantity and may be 
investigated instead. In that case the corresponding local 
density operators read 

B B 
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Thus, according to the scheme in Sec. |TT1 the associated 
current is given by 



1.5 



(37) 



and [J, Ho] = 0. Particularly, (J^) = as well 

as = N{ABl + Jl + ABl)/lQ. 

Because of the commutation of the operators J and 
Hq a second order prediction may be formulated for the 
limit of weak 'interactions', i.e., small z-components of 
the magnetic field. In fact, direct numerics for ~ 16 
by the use of ED already indicate a well-behaved second 
order decay rate R2{t), i.e., i?2(i) appears to take on a 
constant value at long time scales, see Fig. [HI Although 
TV ~ 16 is still far away from the thermodynamic limit, 
the convergence with A'' in Fig. |9] seems to be at least 
rather indicative for a constant value. However, since at 
short time scales i?2(i) shows a non-trivial dependence on 
time, such a dependence certainly occurs also for higher 
order decay rates. In particular higher order decay rates 
may not develop towards constant values at long time 
scales. On that account a second order prediction has to 
be considered carefully here, see Sec. IIIII 
Nevertheless, the convergence in Fig. |9] is sufficient for 
strong interactions^^., analogously to the previous case of 
an alternating magnetic field. The resulting second order 
prediction for is shown in Fig. [TU] for equally large 
components of the magnetic field, i.e., B^/B^ ~ 1. The 
use of the fourth order approximation again allows to 
correctly describe 'D{t) up to the (first) zero-crossing of 
the underlying current autocorrelation function^. But, 
in contrast to the decrease in Fig. |8l a renewed increase 
of 'D{t) emerges after this zero-crossing, resulting from 
partial revivals of the current autocorrelation function, 
e.g., due to a spectrum which gradually becomes closer 
to the equidistant levels of the pure Ising model (with a 
magnetic field in z-direction) . This renewed increase of 
I?(t) remarkably turns out to be captured by the mere 
second order prediction. Therefore the present example 
clearly illustrates that fourth order corrections improve 
usually the description at short time scales but do not 
yield necessarily to a better description at larger time 
scales, e.g., after a potential zero-crossing of the current 
autocorrelation function. 
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FIG. 9. (color online) The second order decay rate R2{t) for 
the energy current in the Ising chain in the presence of a tilted 
magnetic field for A'^ = 10 and 16. Parameters: Bx/Jcc = 
0.85. 
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FIG. 10. (color online) The time-dependent diffusion constant 
T>{t) for energy transport in the Ising chain in the presence of 
a tilted magnetic field for (a) Bx/B^ = 1.0 (a = 45°) and (b) 
-Ba;/Bz = 0.85 {a 40°). Parameters: B^/Jcc = 0.85. The 
numerical results from ED (black, solid curves) are shown for 
N = 10, 14, and 16. The theoretical predictions according to 
TCL2 (red, dashed curves) are also based on numerical data 
for iV = 16 in Fig. H 



VII. SUMMARY AND CONCLUSION 

The present paper has studied the decay of current 
autocorrelation functions for quantum systems featuring 
strong 'interactions'. In this study the term interaction 
has referred to that part of the Hamiltonian causing the 
(major) decay of the current. To this end an appropriate 
perturbation theory in the interaction strength has been 
introduced at first, namely, by an application of the TCL 
projection operator technique. For the addressed case of 
strong interactions the quality of a truncation to lowest 



order has been demonstrated to depend on the form of 
the interaction and not on its strength as such. By the 
use of the introduced perturbation theory the diffusion 
coefficient has been evaluated afterwards for a variety of 
transport quantities in concrete quantum systems. This 
evaluation has been started for excitation transport in 
the modular quantum system and has been continued 
for spin and energy transport in several spin chains. For 
all examples the lowest order prediction for the diffusion 
constant has well agreed with the numerical results from 
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ED, even in the case of strong interactions. Remarkably, 
higher order corrections have played a minor role. 

The investigation has focused on high temperatures 
and one-dimensional quantum systems so far. Both have 
been chosen here in order to allow for a comparison with 
the numerical results from ED, being more or less free 
of finite size effects for this particular choice. However, 
the introduced perturbation theory is not restricted to 
one dimension. But this perturbation theory probably is 
restricted to high temperatures, at least in the addressed 
case of strong interactions. Low temperatures confine the 
perturbation theory in the form at hand to the case of 
weak interactions: Only for that case an approximation 
of the statistical operator on the basis of the uncoupled 
system is expected to be reliable at all. Nevertheless, in 
the context of weak interactions, higher order corrections 



probably play a more important role, similarly to high 
temperatures. In any case further estimations for higher 
order contributions certainly are desirable, not only for 
the projection onto currents but also for the alternative 
projection onto densities. These projections have already 
been applied successfully for one-particle models^SiS^. 
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